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DEAD ZONE IN THE POLAR-CAP ACCELERATOR OF PULSARS 
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ABSTRACT 

We study plasma flows above pulsar polar caps using time-dependent simulations of plasma particles 
in the self-consistent electric field. The flow behavior is controlled by the dimensionless parameter 
a = j /cpcj where j is the electric current density and /5qj is the Goldreich- Julian charge density. 
The region of the polar cap where < a < 1 is a "dead zone" — in this zone particle acceleration 
is inefficient and pair creation is not expected even for young, rapidly rotating pulsars. Pulsars with 
polar caps near the rotation axis are predicted to have a hollow-cone structure of radio emission, as 
the dead zone occupies the central part of the polar cap. Our results apply to charge-separated flows 
of electrons (j < 0) or ions (j > 0). In the latter case, we consider the possibility of a mixed flow 
consisting of different ion species, and observe the development of two-stream instability. The dead 
zone at the polar cap is essential for the development of an outer gap near the null surface pc.T = 0. 
Subject headings: plasmas — stars: magnetic fields, neutron 



1. INTRODUCTION 

Magnetic field lines that pass through the light cylinder 
of a rotating neutron star are twisted and carry electric 
currents js = (c/47r)VxB. These currents are sustained 
by electric field En induced along the magnetic field B, 
and ohmic dissipation E\\j feeds the observed pulsar ac- 
tivity. Voltage associated with E\\ controls the energies 
of accelerated particles, creation of secondary electron- 
positron pairs, and emission of radio waves. The accel- 
erating voltage has been discussed in a number of works 
on pulsars beginning from early papers in the 1970s 
(Sturrock 1971; Ruderman & Sutherland 1975; Arons 
& Scharlemann 1979). 

The key dimensionless parameter of the polar-cap ac- 
celerator is 

3B n x 

a= , (1) 

CPGJ 

where pgj = — O • B/2-7TC is the local corotation charge 
density of the magnetosphere (Goldreich & Julian 1969). 
For a special value of a — a (close to unity) a steady 
state was found for the polar-cap flow with significant 
particle acceleration (e.g. Arons & Scharlemann 1979; 
Muslimov & Tsygan 1992). However, a is not, in gen- 
eral, expected to take this special value (e.g. Kennel et al. 
1979). Global solutions for approximately force- free pul- 
sar magnetospheres give a that significantly varies across 
the polar cap (Timokhin 2006). In general, a can take 
any value from — oo to +oo, depending on the polar cap 
distance from the rotation axis and the location inside 
the polar-cap region. 

The character of the polar-cap accelerator strongly de- 
pends on a (Mestel et al. 1985; Beloborodov 2008, here- 
after B08). The steady solution with a = cto » 1 is a 
separatrix between two opposite regimes of efficient and 
inefficient acceleration^] In particular, if < a < 1, Eu 
is quickly screened in the charge-separated plasma flow- 

1 Hereafter we will refer to this separatrix as a = 1, neglecting 
the deviation of cto from unity. The precise ciq is controlled by the 
curvature of magnetic field lines and the general relativistic effects 



ing from the polar-cap surface. The electric field satisfies 
Maxwell equations that read (in the co-rotating frame of 
the star, see e.g. Fawley et al. 1977; Levinson et al. 
2005), 

V • E = 4tt(p - poj), (2) 



<9E 



= Mj B -j). 



(3) 



If < a < 1, there exists a velocity v = ac that allows 
the charge-separated flow j = pv to simultaneously sat- 
isfy p = pgj and j = js- If the flow started from the 
conducting boundary (which has E — 0) with v = ac, 
no electric field would be generated (then V • E = and 
d~E/dt — 0). The actual boundary has v ^ ac, as charges 
are lifted from the polar-cap surface with a small initial 
v, comparable to the thermal velocity in the surface ma- 
terial. The deviation of v from ac implies p ^ pej or 
j / jg, which generates electric field. B08 argued that 
Equations ([2]) and ([3|) with < a < 1 always drive the 
flow toward v — ac, like a pendulum is driven by gravity 
toward its equilibrium position. The resulting oscilla- 
tions occur in space or time, according to Equations @ 
or ([3]), respectively. For example, the steady-state solu- 
tion for a cold flow exhibits oscillations in space (Mestel 
et al. 1985; B08). The oscillatory behavior of the flow 
with < a < 1 is, in essence, Langmuir oscillations; 
they are generated near the boundary where the flow is 
injected with v < ac and accelerated toward v = ac. 

In this paper, we investigate the accelerator with < 
a < 1 in more detail. In Section 2, we write down 
the steady-state solution for the charge-separated flow, 
generalized to non-zero temperature of the polar-cap. 
We argue that the flow is unstable to small perturba- 
tions and can develop into a complicated time-dependent 
state with a broad momentum distribution. To explore 
the behavior of the flow, we perform fully kinetic time- 
dependent simulations. The method of simulations is de- 

(Muslimov & Tsygan 1992); its exact value is close to unity and is 
not essential for the rest of the paper. 
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scribed in Section 3, and the results are presented in Sec- 
tions 4 and 5. Our simulations confirm the predicted tur- 
bulent Langmuir oscillations with a small voltage. Parti- 
cle acceleration in the flow with < a < 1 is insufficient 
to ignite pair creation. Implications of this "dead zone" 
for radio emission and outer gaps in pulsars are discussed 
in Section 6. 

2. STEADY-STATE SOLUTION FOR A 
CHARGE-SEPARATED FLOW 

2.1. Basic equations 

It is natural first to attempt to construct a simple 
model assuming that the polar-cap flow is steady in the 
(rotating) frame of the neutron star. Given the steady 
magnetic field in this frame, and the steady boundary 
conditions at the stellar surface — excellent static con- 
ductor that can supply charges with a given temperature, 
— one could expect a steady state to be established un- 
less the flow is prone to an instability. 

Consider a charge-separated flow from the polar cap 
that carries electric current js along magnetic field B. 
In a steady state j = jb (Equation [3]). For simplicity, let 
us assume that B is approximately perpendicular to the 
polar cap and let z measure the altitude above the stellar 
surface. A particle of mass m and charge e that starts 
with a Lorentz factor 70 ~ 1 at z — will accelerate as 
it moves along the magnetic field line, 



l{ z ) = 7o +a(z), 



e($-$ ) 



mc 



(4) 



where $ is the electric potential and £ji = —d&Jdz. 
Gravitational acceleration (and centrifugal acceleration 
in the rotating frame) is neglected compared to the elec- 
tric acceleration. 
The electric potential satisfies Poisson equation, 



lb 2 



An(p - p GJ ), 



(5) 



where we assumed that the potential varies along z much 
faster than it does in the transverse directions, i.e. the 
acceleration length lu is much smaller than the charac- 
teristic transverse scale of the problem l±, which may 
be associated with the size of the polar cap. This con- 
dition is satisfied for the flows considered below0 The 
term — /?gj may be viewed as a fixed background charge 
density. The charge density of the flow itself is given by 



p(z) = 3b J 



w(7o) cfyp 
"(70,2) 



(6) 



Here u(7o, z) is the velocity of particles that start at z — 
with initial Lorentz factor 70; note that v 2 /c 2 = 1 — j~ 2 
where 7(2;) is given by Equation (j4j. Function 1^(70) 
describes the probability distribution of 70 . The width 
of this distribution is controlled by the temperature of 
polar cap T. For example, w = $(70 — 1) describes a cold 
polar cap (T = 0) where all particles have 70 = 1. 

2 Alternatively, the additional term V^"!? could be moved to 
the right-hand side of Equation 10 and included in the effective 
PGJ- 



We multiply both sides of Equation §5$ by da/dz 
-(e/mc 2 )d$/dz, substitute Equation ©, and find 



'da V 
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dp . da 

"Two, z ) w \1o) "7o - ~r Pgj 
dz dz 

(7) 



On the right-hand side, we used da/dz = —d'y/dz (Equa- 
tion [4]) and d~f/v = dp/c. Integration of Equation (J7J in 
z gives 



A? / da\' 



2 \dz J 
where 



- = 



[p(70,«) -po] w(7o)d7o 



a{z) 



, (8) 



P 2 (7o, z) = 7 2 - 1 = [70 + a(z)Y 



(9) 

In Equation ((5} we used a(0) = and the boundary 
condition da/dz(0) — (the stellar surface is modeled 
as a perfect conductor that can freely emit charges with 
£7||(0) = 0). We also used Jb(z) rs const and pgj{ z ) ~ 
const, as jb and pqj do not significantly vary on the 
characteristic acceleration length X p , which is defined by 



^2 mc 
p AirejB 



(10) 



This length may be thought of as the plasma skin depth; 
it is related to the plasma frequency Li p , 



An 



Aime 2 



(11) 



where n — js /ec is the characteristic plasma density. 

A quick estimate for js and X p in pulsars may be ob- 
tained from the following consideration. The magnetic 
flux through the polar cap 'J equals the flux through the 
light cylinder _Rlc = c/fi. The bundle of open field lines 
is strongly twisted at the light cylinder (toroidal compo- 
nent comparable to poloidal), and hence it carries electric 
current / ~ c^/27ri?LC: according to Stokes theorem. 
The current density near the star satisfies Jb/B ~ I/*f> 
(which follows from the fact that j flows along B); this 
yields 

VLB , x 

is - ^- (12) 

Then the plasma skin depth in the polar-cap accelerator 
may be expressed as 



Ar, 



eB 
mc 



(13) 



The scale X p is much smaller than the typical size of the 
polar cap r pc ~ (R^Q/c) 1 / 2 , where i?NS ~ 10 6 cm is 
the radius of the neutron star. 

2.2. Cold and warm solutions 

Once the injection distribution 10(70) is specified, it 
is straightforward to numerically integrate Equation (|8]l 
and find a(z). In our sample models we chose 10(70) = 
(fcT)- 1 exp[-(7 - l)/kT] with kT/mc 2 = (cold) and 
0.03 (warm); the average injection momentum po in the 
warm model equals 0.22mc. Figure [T] shows $(z) for the 
cold and warm solutions. 
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Fig. 1. — Steady-state solution for the charge-separated polar- 
cap flow with a = 0.8. Two cases are shown: cold polar cap 
T = (solid curve) and hot polar cap kT/mc 2 = 0.03 (dashed 
curve), which corresponds to average injection momentum 0.22mc. 
Dotted curve shows the solution for a flow where all particles are 
injected with the same po = 0.22. 



Figure [T] also shows a third model where all particles 
injected at the polar cap have po = 0.22, i.e. 10(70) is a 
delta-function. In this model, Equation © simplifies to 



\l ( da\ 2 a(z) 

T =P(lo,z) — po H : 

dz I a 



(14) 



[same as Equation (3) in B08]. This flow is everywhere 
cold, i.e. its momentum distribution is described by 
f(p') = n-S\p' — p(z)]. As one can see in Figure [TJ the 
cold model with po 7^ provides an excellent approxima- 
tion to the exact warm model that has the same average 
value of po ■ 

The cold flow solution was discussed in earlier works 
(Mestel et al. 1985; B08). For 1 - a < 1, the oscillation 
period is approximately given by (B08) 



zo 



,3/2. 



Ar, 



1 



1 - a < 1. 



(15) 



The precise period is obtained by numerical integration; 
e.g. z = 11.0 X p for a — 0.8. The momentum of the 
steady cold flowp(z) oscillates between the injection mo- 
mentum po <C 1 and a maximum value p max - The min- 
ima and maxima are where da/dz = 0, and from Equa- 
tion (| 14|) one finds 



Pu 



2a 7o - (1 + a 2 )p 



(16) 



The above equations assumed a(z) = const. In real 
pulsars, a varies due the field-line curvature and gen- 
eral relativistic effects (Muslimov & Tsygan 1992). The 
length-scale of this variation (da/dz) -1 is typically larger 
or comparable to the stellar radius, which exceeds X p 



by several orders of magnitude. When a varies with z, 
the analytical integration of the dynamic equation is not 
possible and one has to solve the two coupled differential 
equations (B08), 



(17) 
(18) 




a(z) 



The solution is similar to the case where a is constant, 
as long as < a < 1. The momentum p(z) quasi- 
periodically passes through maxima and minima. The 
only difference is that the period zq and p max now grad- 
ually change with z (see Figure 1 in B08). 

2.3. Stability of the flow 

Although the cold flow solution with po = 0.22 repro- 
duces very well the electric potential &(z) of the exact 
warm solution with the same average po, the warm and 
cold flows are qualitatively different. Their different mo- 
mentum distribution functions f{p,z) leads to a qualita- 
tively different response to small perturbations. 

Consider first the cold-flow solution shown by the blue 
dotted curve in Figure [T] Since all particles are injected 
with the same momentum po = 0.22, all of them follow 
a single trajectory in the phase space (z,p). They peri- 
odically reach the minimum momentum po at z£ = kzo 
(k = 0, 1, ...) where potential $ reaches maximum. There 
are no particles with momenta p « 0, so a small perturba- 
tion cannot force any particles to reverse their direction 
of motion, and hence the perturbation will be advected 
along the flow. This flow is expected to be stable. 

In contrast, the warm flow (dashed curve in Figure [T]) 
has a broad distribution of po that extends from po = 0. 
At each peak of the electric potential (at z = z%) there 
is a population of particles with nearly zero velocities. 
Consider a perturbation at z w z\. For example, sup- 
pose a small bunch A of particles with momenta in a 
range (pi , pi + Ap) are slightly pushed forward while the 
rest of particles are unperturbed. This perturbation im- 
plies a local increase in electric current j > jg and hence 
dEu/dt < (Equation [3]), generating negative electric 
field SEu at z « z\ that tends to restore the condition 
j = js- In contrast to the initial perturbation, the in- 
duced 5E\\ affects all local particles, regardless of their 
momenta, not just bunch A. This has two implications: 
(1) The induced Eu < will easily and quickly reduce j 
back to js but will be unable to decelerate bunch A to 
the momentum it would have in the steady state flow — 
bunch A will continue to move to z > z£ with a larger 
momentum. (2) SEu < will give very slow particles 
p « negative velocities, creating a new bunch B that 
slides backward down the potential hill. Bunch B creates 
j < ]b at z < z* k , and the system reacts there by induc- 
ing a small 5E\\ > 0, which accelerates all local particles, 
regardless their momenta, not just bunch B. As a result, 
j quickly recovers to jb , however, bunch B is not stopped 
from moving backward and away from z — z\. 

One concludes that the perturbation creates a per- 
manent damage to the steady state that broadens the 
momentum distribution by creating backflowing parti- 
cles. This perturbation is not advected away along the 
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flow, and can develop further. The backflowing particles 
turn out to be trapped between two peaks of the elec- 
trostatic potential. Further development can be studied 
with kinetic time-dependent simulations; it eventually 
completely destroys the steady state solution. 

3. NUMERICAL SETUP 

Our numerical method is similar to that used by Be- 
loborodov & Thompson (2007, hereafter BT07). The 
plasma is modeled as a large number N ~ 10 6 of indi- 
vidual particles that flow along the magnetic field lines. 
We assume that the magnetic field is fixed in the co- 
rotating frame of the star; thus js and pqj are fixed. 
Then the problem becomes essentially one-dimensional, 
as discussed in detail in BT07. In the present paper, we 
consider only charge-separated flows, with no pair cre- 
ation. Three other differences from the magnetar simu- 
lation in BT07 are as follows: (1) The magnetar problem 
had a 3> 1 (pgj wa s negligible compared with j'b/c); 
in contrast, pqj is crucial for polar-cap flows considered 
here. (2) The presence of gravity was essential for the 
closed-field circuit considered in BT07, where the global 
plasma flow was studied on a scale comparable to the 
radius of the star; in the problem considered here the 
electric fields are screened on a much smaller scale ~ \ p 
and the gravitational acceleration plays no role. (3) The 
flow behavior on the small scales z <C -Rns m &y be stud- 
ied using a smal computational box H <C -Rns with an 
open outer boundary (see below). 

In the absence of pair creation, the flow is composed of 
particles lifted from the surface. In most simulations pre- 
sented below we assume that all particles have the same 
mass m and charge e. The particle motion is described 
by the equation, 



tion, 



dpi _ eE\\(zj) 
dt 



mc 



i = l,..,N, 



(19) 



where pi is the momentum of the i-th particle in units 
of mc, and En(zi) is the self-consistent electric field at 
the particle location Z{. The field is found by integrating 
Gauss law (Equation [5]) along the magnetic field line, 



E\\(zi) = 47r [eN(zi) - pojZi 



(20) 



Here N(zi) is the column density of particles between 
z = and z = z^ and we used the boundary condition 
En(0) — 0, as the material below the stellar surface is 
assumed to be a very good conductor that can emit free 
charges. Divergence of the perpendicular component of 
electric field Ej_ is neglected in Equation ((20)) (see BT07 
for discussion of this approximation). The approxima- 
tion |Vj_ • Ej_| <C \dE\\/dz\ is valid if the characteristic 
scale of the flow acceleration zq is smaller than the trans- 
verse scale l±, which is limited by the polar-cap size r pc ; 
the condition zq <C r pc is satisfied in the dead-zone mod- 
els presented below. We also assume that pcj is approx- 
imately constant on scale Zq. Equations (fT9|) and (|20|) 
in essence describe a relativistic, time-dependent diode 
problem with an additional fixed background charge den- 
sity -pcj- 

As we track the motion of all particles individually, 
the continuity equation is automatically satisfied; for a 
charge-separated flow it is equivalent to charge conserva- 



^ + ^=0. 
dt dz 



(21) 



Equation © follows from Equations @ and (|2ip . so we 
will not need Equation @. Instead, the parameter js en- 
ters the problem as a boundary condition. The magnetic 
field lines are frozen in the excellent conductor below the 
stellar surface, which sustains j(0) = ]b- This condition 
is enforced in the simulation by injecting the charges in 
the computational box at z = with the fixed rate jb 
(BT07). 

The electric current js is enforced at one boundary 
z = 0. Since the computational box has a finite size H, 
we also have to choose a boundary condition at z = H 
and the value of H. In all sample models shown in this 
paper we use the simplest boundary condition: particles 
moving out of the box are lost and no particles enter the 
box at z = H. This condition may be refined by al- 
lowing a small inflow of returning particles at the outer 
boundary. We ran test simulations that show that the 
refinements are not important as long as the boundary is 
sufficiently far, so that H is much larger than the char- 
acteristic scale of the flow acceleration. 

In the one-dimensional model, the transverse gradients 
are neglected and the flow effectively has a slab geometry. 
Then it is sufficient to follow particles flowing through 
a small area A of the slab. This allows one to chose 
a reasonable number of particles in the computational 
box, N ~ AHn, e.g. N ~ 10 6 , so that their dynamics 
can be followed in a reasonable computational time. On 
the other hand, N should be large enough so that the 
plasma scale X p contains many particles N p — AX p n. 

In summary, we choose N and H so that 



H 



> 1, 



N.„ 



A 



-iV> 1. 
H 



(22) 



In this limit, the results are expected to be independent 
of the choice of N and H (we verified this by varying the 
two parameters). For most of our simulations H — 100 A p 
and N ~ 10 6 . Another requirement is a small time step 
of the simulation, At <C w^ 1 , so that plasma oscillations 
are well resolved. 



4. RESULTS 

4.1. Steady state and stability tests 

In our simulations and in reality the plasma above pul- 
sar polar caps is collisionless. In the absence of pair cre- 
ation it must satisfy the Vlasov equation, 



OF 
~dt 



dp 



WF + ^.V P F = 0, 



(23) 



where F(t,z,p) is the particle distribution function in 
phase space. The electric current is j(t,z) = pv where 
v(t, z) is the average velocity of the particles. As a first 
simple test, consider a uniform flow with p(z) — pgj, 
v(z) — ac, and Eu (z) = 0. It is easy to see from Equa- 
tions (|19|). (|20|) and ([23)) that the flow must remain in 
this state. This behavior is reproduced by our simula- 
tions. The steady uniform flow can have any momentum 
distribution F(p) as long as v = ac. Note that it re- 
quires a continual injection of particles at z = with the 
average velocity v = ac (which also requires < a < 1). 
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As a second test, consider a "cold" flow where all par- 
ticles move with momentum p(z), with zero momentum 
dispersion. Suppose the flow is injected at z — with 
velocity vq < ac. Then E\\ must be generated, accel- 
erating the flow. In a steady state, the solution for the 
cold flow must have the form, F(z,p') = n(z) S[p' — p(z)], 
where p(z) and n{z) can be described analytically. We 
first test the special case a — 1 (Michel 1974). The flow 
is accelerated by the self-consistent E\\{z), and p exceeds 
unity at z ~ X p . At heights z ;g> A p , velocity approaches 
c, charge density of the flow p = js/v approaches pa, 
and electric field E\\ asymptotes to a constant value, 



%-KmcjB 



(7o - Pa) 



1/2 



[l + Oip- 1 )], (24) 



where 70 = (1 — v 2 /c 2 ) x l 2 and p = 7o/?o- Then the 
flow momentum keeps growing linearly with z, 



V{z) = [2(7o 



Pq)\ 



2»A„ 



(25) 



This analytical solution is reproduced by our simulations. 
After an initial relaxation period (comparable to the light 
crossing time of the computational box) the system for- 
got initial conditions and relaxed to the steady state 
shown in Figure [2] (in this example, vo = 1/6). The 
charge density of the flow is large near the polar cap 
surface and asymptotes to pgj at z ^ A p , as expected. 

Then we studied cold flows with < a < 1 with a 
fixed injection velocity /3q. We chose in our sample nu- 
merical model a — 0.8 and /3o = 0.2. The computational 
box was initially empty; the plasma injected at z — 
filled the box on the dynamical timescale ~ H/c and es- 
tablished a steady state shown in Figure El The steady 
state is in perfect agreement with the analytical model of 
Section 2. The charge density p{z) has spikes at z — kzo 
(k = 0, 1, ...) where the flow has the minimum velocity (3$; 
the height of each spike is p max = js/Po = (a//3o)PGJ- 
The charge spikes are associated with maxima of the elec- 
tric potential (Figure|3b). The oscillating momentum has 
maxima p max = 3.6, in excellent agreement with Equa- 
tion P^|) . The period of oscillation is z 
as found using the method of Section 2. 

As anticipated in Section 2.3, we find that the steady 
state becomes unstable if we reduce /3q to zero. Then 
any small perturbation (e.g. due to numerical error) 
completely destroys the steady state; instead, a time- 
dependent state forms, with a broadened momentum dis- 
tribution function. A steady flow with a finite f3 ^ 
can also be destroyed, although in this case a finite, 
sufficiently large perturbation is required. In fact, this 
case provides a better setup for a numerical analysis of 
the instability, as we can control the form of the initial 
perturbation and then observe how it destroys the flow 
that was stable before the perturbation was applied. We 
made such an experiment with the flow with a — 0.8 and 
(3q = 0.2. We applied a perturbation that was localized in 
space and time — a small "kick" Sp was given to all par- 
ticles located in a small region Sz — A p /2; in this experi- 
ment dp had a Gaussian distribution with the mean value 
and dispersion equal to 0.02. We observed the following 
evolution. As the localized perturbation moved along 
with the background flow, it was greatly amplified when 
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Fig. 2. — Test run for a cold-flow model with a = 1 and vo = c/6. 
The flow relaxed to a steady state in the entire box H = 10 2 A P on 
the light-crossing timescale, H/c; the state of the system is shown 
at t = lOH/c. (a) Flow momentum per particle p(z) in units of 
mc. (b) Charge density p(z) in units of Pgj- 



it reached the potential maximum (which corresponds to 
the minimum pq « 0.2 of the steady-state solution, see 
Figure [3]) , and some particles acquired a negative mo- 
mentum, i.e. reversed their direction of motion. Most 
of the reversed particles became trapped between two 
potential maxima, and some of them were able to pene- 
trate even further back, beyond the preceding potential 
peak. The perturbation further spread in the phase space 
and the damage to the initial steady-state solution was 
further amplified with time, in particular near the poten- 
tial maxima. Eventually, the entire flow became strongly 
time-dependent and the regular periodic structure of po- 
tential peaks disappeared. 

The amplification of small (linear) perturbations at the 
potential maximum can be understood as follows. Con- 
sider a particle whose Lorentz factor differs from that of 
the background cold flow by a small Sj. As the parti- 
cle moves along with the flow, its deviation ^7 remains 
constant, because it travels in the same electrostatic po- 
tential of the background flow (cf. Equation Using 
the relation d-y/dp = (3, we find the perturbation of mo- 
mentum bp that corresponds to £7, 



Sp - 



57 
J 



tx F 



(26) 



It grows as the particle (and the background flow) de- 
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Fig. 3. — Cold flow with a = 0.8 and /3 = 0.2 at time t = 
lA5H/c (a) Momentum p (in units of me) . Red line shows the 
maximum value predicted by Equation JTBJ. (b) Charge density, 
(c) Electrostatic potential. Red line shows the minimum value 
predicted by the analytical model of Section 2. 



celerates near the potential maximum; the correspond- 
ing amplification factor fi^ 1 is particularly large if /?o is 
small. 

The generation of backflowing particles at the potential 
peaks z£ plays a key role in disrupting the steady state. 
We also observe that, in a flow with a finite minimum 
velocity j3o > 0, only a finite, sufficiently large pertur- 
bation can destroy the steady state. The perturbation 
would need to steal from particles energy 70 — 1 « /3q/2 
so that they can be reflected by the potential hill. The 
energy gap 70 — 1 stabilizes the flow against infinitesimal 
perturbations, and only a sufficiently strong kick disrupts 
the flow. 

The trapped/backflowing particles have a deteriorating 
effect on the steady state because they are not advected 
away with the flow and instead repeatedly approach the 
same potential peaks, amplifying the perturbations. In 
addition, one can view the trapped particles as extra 
charge that distorts the electric field. Let iVt rap be the 
number of particles trapped between two potential peaks 



z^_ l and z\\ they create electric field E' = AireN\ 



trap 



at 



z > Z%. 

static potential $' 



k . The corresponding distortion of the electro- 



-E'z grows linearly with z and 



becomes significant at sufficiently large z even if A^trap 
is small. The distance z required to produce e^' ~ mc 2 
is z ~ (Np/Nt Tap )\p. This behavior is qualitatively con- 
firmed by our numerical experiments with larger simu- 
lation boxes H — the flow was found to become more 
unstable with increasing H. 

4.2. Time- dependent state with warm particle injection 

In a more realistic model, particles are lifted from the 
polar cap with a thermal velocity dispersion Avq ~ vq. 
The flow still starts with a small velocity v c and 
hence with a large charge density p 3> pcj, which self- 
consistently generates the accelerating electric field. The 
basic acceleration mechanism is the same as for the cold 
flow shown in Figures [5] and [3] However, there is a new 
feature: particles with different initial velocities behave 
differently in the collective electric potential, and the 
charge density p(z) is changed from the cold-flow solu- 
tion, even though Avq <C c. Some particles have v w 
and can reverse their motion in the regions of growing po- 
tential (£j| < 0), which greatly complicates the behavior 
of the distribution function F{z,p). 

In our simulations, we modeled the warm injection by a 
one-dimensional Maxwell distribution, which is a simple 
Gaussian with dispersion Avo equal to the mean value 
%; we chose Co = 0.2c. As initial conditions we took the 
steady-state solution (Section 2). The main parameter 
of the flow is a, and we performed simulations for several 
values of a in the range < a < 1. 

As expected, the steady state was quickly destroyed 
and the flow kept oscillating in space and time. The ba- 
sic parameters of the flow remained, however, similar to 
the steady cold model. The average charge density (av- 
eraged over oscillations) is nearly equal to pcj and the 
average velocity v is nearly equal to ac, so that the con- 
dition j — jb is satisfied. Figure 0] shows the evolution 
of the hydrodynamic velocity v(t) measured at a fixed 
location z\ (we chose z\ = 50 X p , in the middle of the 
computational box; v was calculated by averaging over 
particles inside a small bin around Zi, of width 2X p ). The 
hydrodynamic velocity v(t) oscillates around ac; these 
oscillations have a relatively small amplitude 5v <C v. 

The moderate value of the hydrodynamic velocity does 
not, in principle, exclude acceleration of a fraction of par- 
ticles to much higher energies. We therefore also stud- 
ied the momentum distribution of particles in the flow. 
Figure [5^ shows a random snapshot of the particle dis- 
tribution in the phase space for the flow with a = 0.8. 
We randomly chose 1000 particles between z = and 
z = 100Ap and the figure shows their locations in the two- 
dimensional phase space (z,p). The simulation demon- 
strates the following: 

(1) There is no high-energy tail in the momentum dis- 
tribution. 

(2) At each z, the momentum distribution has a pro- 
nounced narrow peak at p pea k- Thus, a large fraction 
of particles form a cold stream; this fraction is approx- 
imately equal to a (see below). The momentum of the 
cold stream p pea k is above (but comparable to) p max pre- 
dicted by the steady-state model. 

(3) There is a low-energy wing in the momentum distri- 
bution which extends to negative momenta. This broad 
component of the particle distribution has a hydrody- 



7 




200 300 400 500 600 700 800 900 1000 
Time tcj p 

Fig. 4. — Evolution of the hydrodynamical velocity v of the flow 
measured in the middle of the computational box. Three models 
are shown: a = 0.95 (purple), 0.8 (blue) and 0.6 (dark green). In 
all three cases, the time-average value of v equals a (red horizontal 
lines). 
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Fig. 6. — Mean expectation and standard deviation for the fluc- 
tuating hydrodynamical momentum of the flow measured at the 
center of the computational box. Four simulations are shown, with 
box sizes H/X p = 50, 100, 150, and 200; all four simulations have 
the same parameter a = 0.8. Red line shows the maximum mo- 
mentum predicted by the steady state model with a = 0.8. 
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Fig. 5. — Snapshot of 1000 randomly chosen particles in phase 
space for the flow with a = 0.8 and /3rj = 0.2. Red dashed line 
shows the maximum momentum p max for the steady cold solution 
with the same a = 0.8 and /3rj = 0.2. (a) Random snapshot for 
the simulation with box size H = 100A P . (b) Another random 
snapshot of a similar simulation with a larger computational box 
H = 200 X p . 



namic velocity close to zero; these particles are "trapped" 
and do not contribute much to the current density; how- 
ever they make a significant contribution to charge den- 
sity. In our sample model, about 20% of particles reside 
in the broad trapped component, and this fact has a 
simple explanation. From the point of view of the cold 
stream dynamics, the broad component provides a back- 
ground that offsets the effect of vacuum charge density 
Pgj by the fraction of 20%. This fraction approximately 
equals to 1 — a, so that the other particles (fraction w a) 
may move in the cold stream with v « c and carry js 
without the mismatch in charge density that would gen- 
erate strong En . In essence, the broad component with 
backflowing particles allows the plasma to self-organize 
so that the cold stream can keep v k, c. This is in con- 
trast to the steady-state solution in Section 2 where all 
particles formed a stream with a positive velocity »/a, 
which must be periodically decelerated and accelerated. 

(4) The cold stream momentum p pca k fluctuates in time 
(the corresponding curve in Figure [3] moves in time). 
However, the qualitative form of the phase-space distri- 
bution remains similar to that in Figure [5j 

As seen in Figure [5^, the flow momentum p pea k de- 
creases near the outer boundary of the computational 
box z = H . This is an artifact of the boundary condi- 
tion (free escape with no backflow) , which suppresses the 
density of backflowing particles near the boundary. As 
a result, a modest negative electric held is induced near 
the boundary, decreasing p pca k so that the flow carries 
the required electric current jb- For comparison, Fig- 
ure EJd shows a random snapshot of a similar model (in 
the same interval < z < 100A P ) that has twice as large 
computational box, H = 200A p . As we increase H, the 
boundary effect moves away to larger z, affecting the flow 
properties only at z H . The flow structure inside the 
box (away from the boundary) does not depend on H. 

To check whether the flow momentum depends on the 
size of the computational box, we ran several simulations 
with the same a — 0.8 and different box sizes H. In 
each simulation, we measured the fluid momentum p in 
the center of the box (using a bin Az = 2X P ) at time 



t = 100o; p 1 . The results are shown in Figure [51 There is 
no systematic variation in p with the box size; the small 
variations 10%) are consistent with the fluctuations 
of p in time for each model(3 

The polar-cap flows in pulsars extend through altitudes 
z much larger than our box size H. The fact that our re- 
sults are independent of H confirm the expected behavior 
— the plasma keeps oscillating and particle acceleration 
is quenched everywhere as long as the flow satisfies the 
condition < a < 1. We observe a quasi- uniform and 
quasi-steady behavior in the computational box (apart 
from the initial acceleration region of length ~ 10A P ). In 
a realistic polar-cap flow, each segment of length ~ 100A P 
should behave like our computational box. 

We also ran simulations with varying a(z). We ran 
models with da/dz ~ 2 x 10~ 4 A~ 1 (realistically, da/dz 

should be even smaller, da/dz Si R^ 1 , where R is the star 
radius). The results are similar to the case of a = const. 
The maximum momentum remains comparable to that 
given by Equation (Ti7if as long as < a < 1. 

5. MIXED ION FLOW AND TWO STREAM 
INSTABILITY 

If jb > (which is equivalent to pgj > for a > 0), 
the charge-separated flow puled out from the polar cap is 
made of ions. Different ion species may end up in such a 
flow, and they will be accelerated to different velocities. 

The mixed ion flow shares many features with the 
identical-particle model studied in the previous sections. 
Steady state solutions can be obtained using the method 
described in Section 2. Ions with different masses and 
charges move with different hydrodynamical momenta 
and co-exist in a common, periodic electrostatic poten- 
tial. This steady solution is prone to kinetic instability 
similar to that described in Sections 2 and 4. There 
is, however, an important new feature: the ion streams 
with different hydrodynamical momenta are prone to 
two-stream instability. 

To study the behavior of the mixed ion flow we slightly 
change the setup of our numerical simulation. Consider, 
e.g., a mixture of protons and helium nuclei (alpha- 
particles). The particle injection at z = now consists of 
two ions species; they have charges e\ and e 2 = 2ei, and 
masses m\ and mi = 4mi. The two species are injected 
with equal rates N\ — N%. Then alpha-particles carry 
electric current j 2 — e2-/V 2 that is two times larger than 
the proton current ji = e\N\. Thus, j2 = (2/3)js and 
ji = (1/3)j'b are maintained at the boundary. 

To define a characteristic plasma skin depth X p we use 
Equation (ITUl) where we replace e, m,jB by (ei,mi,j'i) 
or, equivalently, by (e2,m2,j2) (note that e-ijilm-i — 
e iJi/ TO i)- The characteristic plasma frequency is defined 
by uj p = c/X p . 

Figure 7 shows a snapshot of the phase-space distribu- 
tion of ions long after the beginning of the simulation. In 
this sample model a = 0.4. The modest value of a (not 

3 Note that the average momentum p does not correspond 
to the average velocity v shown in Figure f4] in the sense that 
p ^ /3(1 — /3 2 ) - 1 / 2 , because of the broad low-energy tail of the 
distribution function. Compared with v, the calculation of p gives 
a higher weight to fast particles due to the additional factor 7 in 
p = 7/3. The average velocity remains close to etc, and the averaged 
momentum is larger than /3(1 — /3 2 ) - 1 / 2 . 
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Fig. 7. — Snapshot of the mixed ion flow with a = 0.4 at time 
t = 10 3 ujp . Top panel shows the phase space distribution, where 
red dots represent protons and blue dots represent helium ions. 
Bottom panel shows the electric field. 



close to unity) implies modest Lorentz factors of parti- 
cles and the fast development of instabilities. The flow 
exhibits the following features: 

(1) One period of the steady state solution is reproduced 
near the injection boundary z = 0. The period zq w 3A p 
agrees with the result from numerical integration of the 
corresponding steady state model. (This feature is stable 
in our sample model because we chose a high injection 
velocity vq ~ 0.4c.) 

(2) At larger z the periodic flow becomes unstable and 
develops into a configuration similar to that in Figure 5, 
except that now we have two cold variable streams. Be- 
sides the cold streams, there is a broad distribution of 
ions with smaller momenta and a negligible hydrody- 
namic velocity. The origin of this broad component of 
trapped particles was discussed in Section 4.2 and plays 
here a similar role — it is self-organized so that the 
streams may move with a relativistic speed without a 
mismatch in charge density. 

(3) Further from the boundary (at z > 20A p ), the two 
streams develop a two-stream instability. The growth 
rate of the instability may be estimated using an ideal- 
ized model of two cold fluids with densities ri\, n 2 and 
velocities t>i,U2- It is straightforward to derive the dis- 
persion relation for Langmuir modes with frequency to 
and wave- vector k (e.g. Melrose 1986); it gives, 



7^(w — kv\) 2 t|( w — fcw 2 ) 2 

where uj\ — Anniel/mi and w| = 47rn2e 2 /m2. Using 
u>\ ?» ll>2 ~ w p and the characteristic values of velocities 
V\ , u 2 from our simulation, we find from Equation (|27[) 
that the most unstable modes have to comparable to uj p 
and their growth rate is T ~ 0.2uj p . In the simulation 
we observe a slightly smaller T. The distance over which 
the Langmuir waves are amplified is roughly 10A p . As a 
result of the instability, the two streams are smeared out 
at large z, in particular the stream of lighter ions. No 
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significant particle acceleration is seen in the simulation. 

6. DISCUSSION 

We have presented detailed one-dimensional time- 
dependent simulations of the plasma flow extracted from 
the polar caps of neutron stars. The simulations provide 
a fully kinetic description of the flow, with self-consistent 
electric field and particle distribution function. In this 
paper, we focused on the regime < a < 1, where a is 
the main parameter of the flow defined by Equation ([T]). 
In agreement with the estimates of B08, we find that the 
particles are accelerated to Lorentz factors, 



and are not capable of igniting pair creation. In this 
sense, flows with < a < 1 are "dead." They are sus- 
tained by a modest voltage, oscillating in space and time. 
Although the simulation is limited to regions close to the 
pulsar surface, the result does not depend on the sim- 
ulation box size, and hence should describe the entire 
polar cap flow, as long as a remains between and 1. 
The parameter a is expected to vary along the magnetic 
field lines, on a scale comparable to the stellar radius; 
we have verified that this variation does not change the 
oscillating behavior of the flow (see also B08). 

The simulations show how a kinetic instability devel- 
ops and disrupts the ideal periodic structure found in the 
analytical models of the dead zone; the instability mech- 
anism is described in Sections 2 and 4. We find that the 
momentum distribution function has two distinct parts 
— a variable "cold stream" and a broad wing at low mo- 
menta, which includes particles flowing backward to the 
polar cap. The fraction of particles in the cold stream is 
approximately equal to a; the remaining fraction 1 — a 
forms the broad component. Even though the flow is 
turbulent, it shows no signs of particle acceleration to 
energies higher than that of the cold stream. 

The value of parameter a depends on the location and 
geometry of the polar cap. A simplest magnetospheric 
configuration is that of a centered dipole. Then the pa- 
rameter a depends on the angle between the magnetic 
and spin axes, £; besides, it varies across the polar cap. 
For nearly aligned rotators (£ » 0), < a < 1 in the 
central part of the polar cap and a < in a ring-shaped 
zone near the edge of the polar cap (Timokhin 2006; Par- 
frey et al. 2012). In this case, the dead zone occupies 
the central part of the polar cap, and e discharge must 
be confined to the ring, matching the phenomenological 
"hollow cone" model of pulsar emission. In contrast, the 
polar cap of an orthogonal rotator (£ w tt/2) has \a\ ^> 1, 
which enables discharge for the entire polar cap. At 
arbitrary misalignment < £ < 7r/2, the values of a 
are provided by global three-dimensional simulations of 
the magnetospheric structure (e.g. Spitkovsky 2006) and 
should play a key role for the geometry of the radio beam. 

We presented our results using plasma skin depth X p 
as a unit of length and particle rest-mass mc 2 as a unit 
of energy. In this form, the results do not depend on 
the charge or mass of the particles extracted from the 
polar cap, as long as the flow is made of identical par- 
ticles. In particular, Equation (|28|) is valid for both 
electron flow (pej < 0) and ion flow (pgj > 0), and 



the phase-space distribution shown in Figure 5 describes 
both cases. Note that the accelerating voltage is pro- 
portional to the particle mass; voltage implied by Equa- 
tion ([28)) is different for ions and electrons by the factor 
of rrii/m e ~ 2 x 10 3 . The relatively high voltage in the 
ion flow, e$ » rriiC 2 (l + a 2 )/(l — a 2 ) is still hardly suf- 
ficient to ignite pair discharge by a seed electron or 
positron. 

The identical-particle model may not hold for an ion 
flow; in this case, new effects may enter the problem. 
Firstly, heavy ions pulled out from the polar cap may 
not be completely ionized and begin to lose electrons as 
they are accelerated and interact with the X-rays above 
the stellar surface; this process effectively creates new 
charges, reminiscent of pair creation (e.g. Jones 2012). 
Secondly, the ion flow may be a mixture of different nu- 
clei which will be accelerated to different Lorentz factors. 
The mixed ion flow is prone to two-stream instability, 
possibly leading to formation of plasma clumps and gen- 
eration of coherent radio emission. In our simulations, we 
observe the expected two-stream instability, however do 
not observe significant structure (clumps) in the turbu- 
lent flow. This may change in three-dimensional simula- 
tions. The frequency of excited waves (comparable to the 
ion plasma frequency) is in the radio band, and coherent 
emission from clumps could create bright coherent emis- 
sion. It remains to be seen whether this mechanism can 
contribute to the pulsar emission. If it does, it would 
create an additional component of the radio pulse. In 
the case of approximately aligned rotator, the additional 
component would be generated in the central region of 
the polar cap, leading to a "hollow cone + core" structure 
of the radio pulse. 

The charge-separated model of the dead zone can be 
modified to include possible backflowing particles from 
distant parts of the open field-line bundle (e.g. from 
a pair-producing outer gap). These particles can con- 
tribute to the current density and also serve as an addi- 
tional background charge density, which may be modeled 
as a contribution to the effective "vacuum" charge den- 
sity — pgj- This would change the effective a (Lyubarsky 
1992; B08), most likely reducing it. 

An outer gap is expected to form in a charge-separated 
flow near the null surface B ■ fl = (Cheng et al. 1986). 
On a given field line, the outer gap will be screened if it is 
loaded by multiple pairs produced by the discharge at 
the polar cap. Thus, the suppression of discharge near 
the field-line footpoint is an essential condition for the 
existence of an outer-gap accelerator. Therefore, one can 
expect an outer gap to form on field lines with footpoints 
in the dead zone. 

We did not simulate in this paper flows with a > 1 or 
a < 0; in these cases particles must be strongly accel- 
erated. This regime leads to e^- discharge that must be 
unsteady, with a significant intermittent backflow (B08). 
A model for oscillating discharge may be constructed in 
hydrodynamical approximation (Levinson et al. 2005), 
however a fully kinetic description is essential, as demon- 
strated by our results for the dead zone, where a signif- 
icant fraction of particles are trapped and form a broad 
wing at low momenta in the distribution function. The 
discharge simulation can be done using our setup of a 
fixed current j = jb at the stellar surface (BT07) and 
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incorporating pair creation. We defer the simulations 
with a > 1 and a < to a future work. 

When this work was completed, the preprint by Timo- 
khin & Arons (2012) came out. They present simulations 
of charge separated flows, using a similar method, and 
the results agree with our results for < a < 1. They 



also consider flows with a < and a > 1, and find a 
strong unsteady e 1 *" discharge confirming the analysis in 
B08. 
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